library(fixest) # to demean variables
library(multiwayvcov) # to cluster SEs
library(sandwich) # to get HC2 SEs
library(stargazer) # to generate tables
library(mvtnorm) # to simulate betas and vcov matrix

rm(list = ls())

# load data
load(file = 'chData.RData')

# Put each dataset in a different object
ch_reform <- ch_list[[1]]
ch_refTerm <- ch_list[[2]]
ch_pre_post <- ch_list[[3]]
ch_refTermMChanges <- ch_list[[4]]
ch_refTermParty <- ch_list[[5]]
ch_refTermMChangesParty <- ch_list[[6]]

# Exclude "Gonzalo Fuenzalida" and 'Javier Macaya' because their district changed after the reform
ch_reform <- subset(ch_reform, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTerm <- subset(ch_refTerm, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_pre_post <- subset(ch_pre_post, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTermMChanges <- subset(ch_refTermMChanges, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTermParty <- subset(ch_refTermParty, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))
ch_refTermMChangesParty <- subset(ch_refTermMChangesParty, !(bill_author %in% c("Gonzalo Fuenzalida", 'Javier Macaya')))


# M1: G*M C1,C2,C3 (all data)
data_m1 <- fixest::demean(women_tb_htun_false_pct + magAnalysis + female + magAnalysisWoman
	~ session, data = ch_refTerm)
data_m1 <- cbind(data_m1, bill_author = ch_refTerm$bill_author)
m1 <- lm(women_tb_htun_false_pct ~ magAnalysis + female + magAnalysisWoman - 1, data = data_m1)
vcov1 <- cluster.vcov(m1, cluster = data_m1$bill_author)
s1  <- sqrt(diag(vcov1))

# M2: G*M C1,C2,C3 (M changes at the middle of C2)
data_m2 <- fixest::demean(women_tb_htun_false_pct + magAnalysis + female + magAnalysisWoman
	~ session + reform, data = ch_refTermMChanges)
data_m2 <- cbind(data_m2, bill_author = ch_refTermMChanges$bill_author)
m2 <- lm(women_tb_htun_false_pct ~ magAnalysis + female + magAnalysisWoman - 1, data = data_m2)
vcov2 <- cluster.vcov(m2, cluster = data_m2$bill_author)
s2  <- sqrt(diag(vcov2))


# M3: G*M C1,C3
data_m3 <- fixest::demean(women_tb_htun_false_pct + magAnalysis + female + magAnalysisWoman
	~ session, data = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022')))
data_m3 <- cbind(data_m3, bill_author = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022'))$bill_author)
m3 <- lm(women_tb_htun_false_pct ~ magAnalysis + female + magAnalysisWoman - 1, data = data_m3)
vcov3 <- cluster.vcov(m3, cluster = data_m3$bill_author)
s3  <- sqrt(diag(vcov3))

# M4: C2
inSample <- names(table(ch_pre_post$bill_author)[table(ch_pre_post$bill_author) == 2])
ch_pre_post <- subset(ch_pre_post, bill_author %in% inSample)
data_m4 <- fixest::demean(women_tb_htun_false_pct + magAnalysis + female + magAnalysisWoman
	~ bill_author, data = ch_pre_post)
data_m4 <- cbind(data_m4, bill_author = ch_pre_post$bill_author)
m4 <- lm(women_tb_htun_false_pct ~ magAnalysis + magAnalysisWoman - 1, data = data_m4)
vcov4 <- cluster.vcov(m4, cluster = data_m4$bill_author)
s4  <- sqrt(diag(vcov4))

# M5: G*M C1,C2,C3 (only repeated legislators)
data_m5 <- fixest::demean(women_tb_htun_false_pct + magAnalysis + magAnalysisWoman
	~ bill_author + session, data = subset(ch_refTerm, allSessions == 1))
data_m5 <- cbind(data_m5, bill_author = subset(ch_refTerm, allSessions == 1)$bill_author)
m5 <- lm(women_tb_htun_false_pct ~ magAnalysis + magAnalysisWoman - 1, data = data_m5)
vcov5 <- cluster.vcov(m5, cluster = data_m5$bill_author)
s5  <- sqrt(diag(vcov5))

# M6: G*M C1, C2, C3 (M changes at the middle of C2) (only repeated legislators)
data_m6 <- fixest::demean(women_tb_htun_false_pct + magAnalysis + magAnalysisWoman
	~ session + reform + bill_author, data = subset(ch_refTermMChanges, allSessions == 1))
data_m6 <- cbind(data_m6, bill_author = subset(ch_refTermMChanges, allSessions == 1)$bill_author)
m6 <- lm(women_tb_htun_false_pct ~ magAnalysis + magAnalysisWoman - 1, data = data_m6)
vcov6 <- cluster.vcov(m6, cluster = data_m6$bill_author)
s6  <- sqrt(diag(vcov6))

# M7: G*M C1,C3 (only repeated legislators)
data_m7 <- fixest::demean(women_tb_htun_false_pct + magAnalysis + magAnalysisWoman
	~ bill_author + session, data = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022') & l10_l18 == 1))
data_m7 <- cbind(data_m7, bill_author = subset(ch_refTerm, session %in% c('2010-2014', '2018-2022') & l10_l18 == 1)$bill_author)
m7 <- lm(women_tb_htun_false_pct ~ magAnalysis + magAnalysisWoman - 1, data = data_m7)
vcov7 <- cluster.vcov(m7, cluster = data_m7$bill_author)
s7  <- sqrt(diag(vcov7))

# M8: C3
data_m8 <- subset(ch_refTerm, session == '2018-2022')
m8 <- lm(women_tb_htun_false_pct ~ magAnalysis + female + magAnalysisWoman, data = data_m8)
vcov8 <- vcovHC(m8, 'HC2')
s8  <- sqrt(diag(vcov8))

# Table E.17
stargazer(m1, m2, m3, m4, m5, m6, m7, m8, 
	se = list(s1, s2, s3, s4, s5, s6, s7, s8), 
	covariate.labels = c('M', 'Woman', 'M x Woman'), 
	add.lines = list(c('FE by Legislative Term', 'Yes', 'Yes', 'Yes', 'No', 'Yes', 'Yes', 'Yes', 'No'),
		c('FE by Reform', 'No', 'Yes', 'No', 'No', 'No', 'Yes', 'No', 'No'),
		c('FE by Legislator', 'No', 'No', 'No', 'Yes', 'Yes', 'Yes', 'Yes', 'No')),
	dep.var.labels = "Bills' Portfolio on Women's Issues (%)",
	omit.stat = c("ser",'f'),
	star.cutoffs = c(0.10, 0.05, 0.01),
	no.space = TRUE,
	single.row = FALSE,	
	notes.append = FALSE,
	title = "Association Between Legislative Portfolio, Gender, and District Magnitude--2014-2021  Cámara de Diputadas y Diputados de Chile",
	label = 't:mag'
)

#################### Substantive Effect ########################
## Scenarios: Women, M = 2
W2 <- data.frame(mag = 2, female = 1, magAnalysisWoman = 2 * 1)
## Scenarios: Women, M = 3 (post-reform)
W3 <- data.frame(mag = 3, female = 1, magAnalysisWoman = 3 * 1)
## Scenarios: Men, M = 2
M2 <- data.frame(mag = 2, female = 0, magAnalysisWoman = 2 * 0)
## Scenarios: Men, M = 3 (post-reform)
M3 <- data.frame(mag = 3, female = 0, magAnalysisWoman = 3 * 0)
# Women and Men Data
Wdata <- rbind(W2, W3)
Mdata <- rbind(M2, M3)

## Simulate betas and vcov
set.seed(123)
sims1 <- rmvnorm(n = 1000, mean = coef(m1), sigma = vcov1)
sims2 <- rmvnorm(n = 1000, mean = coef(m2), sigma = vcov2)
sims3 <- rmvnorm(n = 1000, mean = coef(m3), sigma = vcov3)
sims4 <- rmvnorm(n = 1000, mean = coef(m4), sigma = vcov4)
sims5 <- rmvnorm(n = 1000, mean = coef(m5), sigma = vcov5)
sims6 <- rmvnorm(n = 1000, mean = coef(m6), sigma = vcov6)
sims7 <- rmvnorm(n = 1000, mean = coef(m7), sigma = vcov7)
sims8 <- rmvnorm(n = 1000, mean = coef(m8), sigma = vcov8)

## Calculate Predicted Values
# Women
pred1_W <- apply((sims1 %*% t(Wdata[2,])) - (sims1 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred2_W <- apply((sims2 %*% t(Wdata[2,])) - (sims2 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred3_W <- apply((sims3 %*% t(Wdata[2,])) - (sims3 %*% t(Wdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred4_W <- apply((sims4 %*% t(Wdata[2, -2])) - (sims4 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) # omit female
pred5_W <- apply((sims5 %*% t(Wdata[2, -2])) - (sims5 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) # omit female
pred6_W <- apply((sims6 %*% t(Wdata[2, -2])) - (sims6 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) # omit female
pred7_W <- apply((sims7 %*% t(Wdata[2, -2])) - (sims7 %*% t(Wdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) # omit female
pred8_W <- apply((sims8 %*% t(cbind(1, Wdata)[2,])) - (sims8 %*% t(cbind(1, Wdata)[1,])), 2, quantile, c(0.975, 0.5, 0.025)) # add intercept

# Men
pred1_M <- apply((sims1 %*% t(Mdata[2,])) - (sims1 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred2_M <- apply((sims2 %*% t(Mdata[2,])) - (sims2 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred3_M <- apply((sims3 %*% t(Mdata[2,])) - (sims3 %*% t(Mdata[1,])), 2, quantile, c(0.975, 0.5, 0.025))
pred4_M <- apply((sims4 %*% t(Mdata[2, -2])) - (sims4 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) # omit female
pred5_M <- apply((sims5 %*% t(Mdata[2, -2])) - (sims5 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) # omit female
pred6_M <- apply((sims6 %*% t(Mdata[2, -2])) - (sims6 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) # omit female
pred7_M <- apply((sims7 %*% t(Mdata[2, -2])) - (sims7 %*% t(Mdata[1, -2])), 2, quantile, c(0.975, 0.5, 0.025)) # omit female
pred8_M <- apply((sims8 %*% t(cbind(1, Mdata)[2,])) - (sims8 %*% t(cbind(1, Mdata)[1,])), 2, quantile, c(0.975, 0.5, 0.025)) # add intercept

# Building blocks of table 2
bill_hat_women <- c(round(pred1_W[2,], 2), round(pred2_W[2,], 2), round(pred3_W[2,], 2),
	round(pred4_W[2,], 2), round(pred5_W[2,], 2), round(pred6_W[2,], 2),
	round(pred7_W[2,], 2), round(pred8_W[2,], 2))
ci_women <- c(paste0('(', round(pred1_W[3,], 2), ', ', round(pred1_W[1,], 2), ')'), 
	paste0('(', round(pred2_W[3,], 2), ', ', round(pred2_W[1,], 2), ')'), 
	paste0('(', round(pred3_W[3,], 2), ', ', round(pred3_W[1,], 2), ')'), 
	paste0('(', round(pred4_W[3,], 2), ', ', round(pred4_W[1,], 2), ')'), 
	paste0('(', round(pred5_W[3,], 2), ', ', round(pred5_W[1,], 2), ')'), 
	paste0('(', round(pred6_W[3,], 2), ', ', round(pred6_W[1,], 2), ')'), 
	paste0('(', round(pred7_W[3,], 2), ', ', round(pred7_W[1,], 2), ')'), 
	paste0('(', round(pred8_W[3,], 2), ', ', round(pred8_W[1,], 2), ')'))


bill_hat_men <- c(round(pred1_M[2,], 2), round(pred2_M[2,], 2), round(pred3_M[2,], 2),
	round(pred4_M[2,], 2), round(pred5_M[2,], 2), round(pred6_M[2,], 2),
	round(pred7_M[2,], 2), round(pred8_M[2,], 2))
ci_men <- c(paste0('(', round(pred1_M[3,], 2), ', ', round(pred1_M[1,], 2), ')'), 
	paste0('(', round(pred2_M[3,], 2), ', ', round(pred2_M[1,], 2), ')'), 
	paste0('(', round(pred3_M[3,], 2), ', ', round(pred3_M[1,], 2), ')'), 
	paste0('(', round(pred4_M[3,], 2), ', ', round(pred4_M[1,], 2), ')'), 
	paste0('(', round(pred5_M[3,], 2), ', ', round(pred5_M[1,], 2), ')'), 
	paste0('(', round(pred6_M[3,], 2), ', ', round(pred6_M[1,], 2), ')'), 
	paste0('(', round(pred7_M[3,], 2), ', ', round(pred7_M[1,], 2), ')'), 
	paste0('(', round(pred8_M[3,], 2), ', ', round(pred8_M[1,], 2), ')'))

tableE18 <- rbind("Women "= bill_hat_women, " " = ci_women, 
	"Men "= bill_hat_men, " " = ci_men)
colnames(tableE18) <- paste0('Model ', 1:8)
stargazer(tableE18,
	title = "Predicted percentage point change in the cosponsorship portfolio on women’s issues when M increases by one")
